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A Computationally Optimal Randomized Proper 
Orthogonal Deeomposition Teehnique 

Dan Yu, Suman Chakravorty 


Abstract —In this paper, we consider the modei reduction 
prohiem of iarge-scaie systems, such as systems obtained through 
the discretization of partiai differentiai equations. We propose a 
computationaiiy optimai randomized proper orthogonai decom¬ 
position (RPOD*) technique to obtain the reduced order modei by 
perturbing the primai and adjoint system using Gaussian white 
noise. We show that the computations required by the RPOD* 
aigorithm is orders of magnitude cheaper when compared to 
the baianced proper orthogonai decomposition (BPOD) aigorithm 
and BPOD output projection aigorithm white the performance 
of the RPOD* aigorithm is much better than BPOD output 
projection aigorithm. It is computationaiiy optimai in the sense 
that a minimai number of snapshots is needed. We aiso reiate the 
RPOD* aigorithm to random projection atgorithms. The method 
is tested on two advection-diffusion probiems. 

Keywords—Model Reduction, Proper Orthogonal Decomposition 
(POD), Randomization Algorithm. 

I. Introduction 

In this paper, we are interested in the model reduction 
of large scale systems such as those governed by partial 
differential equations (PDE). The dimension of the system 
is large due to the discretization of the PDEs. Eor instance, 
consider the atmospheric dispersion of an air pollutant |[T|. The 
emission of the contaminants on the ground level is shown in 
Fig.IU with four point sources labeled from S1 to S4. 

This is a three dimensional problem, and after discretizing 
the PDE, the dimension of the system is 10®. Therefore, we are 
interested in constructing a reduced order model (ROM) that 
can capture the input/output characteristics of the large model 
such that this ROM can then be used by a hltering algorithm 
for updating the states of the field, such as the Kalman hlter. 
Also, the actuators and sensors could be placed anywhere in 
this held, which leads to a model reduction problem of a large- 
scale system with a large number of inputs/outputs. There 
are two popular contemporary model reduction techniques that 
have been studied in the past few decades; Principal component 
analysis (PCA) and randomization algorithms. 

Among the PCA algorithms. Balanced Proper Orthogonal 
Decomposition (BPOD) ||2l, |0 based on the balanced trun¬ 
cation a and the snapshot Proper Orthogonal Decomposi¬ 
tion (POD) technique 0 has been widely used. Balancing 
transformations are constructed using the impulse responses 
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Fig. 1. Air pollutant problem 


of both the primal and adjoint system, and hence, the most 
controllable and observable modes can be kept in the ROM. 
In 1978, Rung 0 presented a new model reduction algorithm 
in conjunction with the singular value decomposition (SVD) 
technique, and the eigensystem realization algorithm (ERA) 
0 was developed based on this technique. The BPOD is 
equivalent to the ERA procedure ||8l, and forms the Hankel 
matrix using the primal and adjoint system simulations as 
opposed to the input-output data as in ERA. More recently, 
there has been work on obtaining information regarding the 
dominant modes of the system based on the snapshot POD, 
followed by an eigendecomposition of an approximating linear 
operator, called the dynamic mode decomposition (DMD) 0, 
cni. In ca, im, an adaptive POD algorithm based on the 
snapshot POD algorithm is introduced to recursively revise the 
locally valid ROMs. When there are new snapshots collected, 
decisions are made if the basis functions need to be updated. 
Both the DMD and adaptive POD algorithm are applicable for 
the nonlinear systems. 

The primary drawback of the BPOD and ERA is that for 
a large scale system, such as that obtained by discretizing a 
PDE, with a large number of inputs/outputs, the computational 
burden incurred is very high. There are two main parts to 
the computation; the hrst is to collect datasets from compu¬ 
tationally expensive primal and adjoint simulation in order to 
generate the Hankel matrix. The second part is to solve the 
SVD problem for the resulting Hankel matrix. 

To reduce the computational cost of BPOD, improved al¬ 
gorithms have been proposed. Eor example, 0 proposed an 
output projection method to address the problem when the 
number of outputs is large. The outputs are projected onto a 
small subspace via an orthogonal projection Ps that minimizes 
the error between the full impulse response and the projected 
impulse response. However, the method cannot make any claim 
regarding the closeness of the solution to one that is obtained 
from the full Hankel matrix, and is still faced with a very 
high computational burden when both the numbers of inputs 
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and outputs are large. There have also been methods proposed 
0 to reduce the number of snapshots, however, the primary 
problem regarding large number of inputs/outputs remains the 
same. 

There are two major classes of randomization algorithms 
used for low-rank matrix approximations and factorizations: 
random sampling algorithms and random projection algo¬ 
rithms. For a large scale matrix H, random sampling algo¬ 
rithms construct a rank k approximation matrix H by choosing 
and rescaling some columns of H according to certain sam¬ 
pling probabilities Ilia, so the error satisfies \\H — iT||_F ^ 
\\H — ||f + IIj’, with high probability, where is 

a best rank k approximation of H, e is a specified tolerance, 
and ||iF||F denotes the Frobenius norm of H. The improved 
algorithm proposed in ifTSl is to sample some columns accord¬ 
ing to leverage scores, where the leverage scores are calculated 
by performing the SVD of H, so that the error satisfies 
\\H — H\\f < (1 + ^)\\H — iF(fc)||F, with high probability. 
A direct application of both algorithms would require the full 
Flankel matrix to be constructed, however, such a construction 
of the Hankel matrix is computationally prohibitive when the 
number of inputs/outputs is large. Further, the leverage scores 
are calculated by performing the SVD of the Hankel matrix, 
which is also computationally prohibitive owing to the size of 
the problem. 

In random projection method ifThl . the large matrix H is 
projected on to an orthonormal basis Q such that the error 
satisfies ||iJ — QQ'H\\ < (1 -f e)\\H — iF(/c)|| with high 
probability, where ||7J|| denotes the spectral 2-norm of H, (.)' 
denotes the conjugate transpose of (.). A Gaussian test matrix 
fl is generated, and the orthonormal basis Q is constructed by 
performing a QR factorization of the matrix product Hfl. The 
bottleneck of this algorithm remains, as above, the construction 
of the full Hankel matrix, which is prohibitively expensive. 

We had introduced an RPOD algorithm in ifTTl that ran¬ 
domly chooses a subset of the input/output trajectories. A sub- 
Hankel matrix is constructed using these sampled trajectories, 
which is then used to form a ROM in the usual BPOD fashion. 
The Markov parameters of the ROM constructed using the 
sub-Hankel matrix were shown to be close to the Markov 
parameters of the full order system with high probability. 
We proved that a lower bound exists for the number of the 
input/output trajectories that need to be sampled. The major 
problem associated with this algorithm is that different choices 
of the sampling algorithms would lead to different lower 
bounds, and the choice of a good sampling algorithm other 
than the uniform distribution is unclear. 

In this paper, we propose the RPOD* algorithm which is 
closely related to the random projection algorithm. In the 
RPOD* algorithm, we perturb the primal and adjoint system 
with Gaussian white noise, and we prove that similar to 
the BPOD algorithm, the controllable and observable modes 
are retained in the ROM. The Markov parameters of the 
ROM constructed using RPOD* are shown to be close to the 
Markov parameters of the full order system, while the error is 
bounded. The main contribution of the RPOD* method is that 
it is computational orders of magnitudes less expensive when 


compared to the BPOD/ERA and randomization algorithms for 
a large-scale system with a large number of inputs/outputs. The 
RPOD* algorithm can be viewed as applying the random pro¬ 
jection on the full Hankel matrix H twice without constructing 
the full Hankel matrix H, i.e., H = 

where Hi,H 2 are two random projection matrices, and Z,X 
are the usual impulse response matrices of the adjoint and 
primal system. However, we actually only generate Zfl 2 and 
Xfli which can be constructed from a single white noise 
perturbed trajectory each of the adjoint and primal system 
respectively, and thus, are orders of magnitude smaller in size 
than the impulse responses Z and X. Thus, the computational 
cost to generate the Hankel matrix and to solve the SVD 
problem is saved by orders of magnitude. We believe that it is 
the most computational efficient POD algorithm. In practice, 
the RPOD* algorithm can be solved in real-time. 

The rest of the paper is organized as follows. In Section 
im the problem is formulated, and in Section [Till we review 
the BPOD algorithm and illustrate in a simplified fashion how 
to relate the BPOD ROM to the controllable and observable 
modes of the system. The RPOD* algorithm is proposed in 
Section HV] and the formal proofs and results are shown. Also, 
we discuss some implementation problems in this section. In 
Section |V] we compare the RPOD* algorithm with BPOD, 
random projection and BPOD output projection algorithm. In 
Section |Vll we provide computational results comparing the 
RPOD* with the BPOD/BPOD output projection for a one 
dimensional heat transfer problem and a three dimensional 
atmospheric dispersion problem. 

H. Problem Formulation 

Consider a stable linear input-output system: 

Xk = Axk-i + Buk, 

yk=Cxk, (1) 

where Xk S , Uk C 3?^, Uk C are the states, inputs, and 
outputs at discrete time instant tk respectively. Assume that 
A, B, C matrices are real. 

The adjoint system is defined as: 

Zk — A Zk — 1 C Vk , Wk — B Zk , (2) 

where Zk G Wk G 3?^ is the state and output of the adjoint 
system at time tk respectively, Ufe G 

Definition 1. The Markov parameters of the system is 
defined as CA^B, i = 1, - ■ ■ ■ 

Assumption 1. Assume that A is diagonalizable and stable. 

Under Assumption 1, let, 

A = VAU', (3) 

where A are the eigenvalues, (U, U) are the corresponding right 
and left eigenvectors. 

Definition 2. A mode {Ai,Ui,Vi) is not controllable if 
C7'i? = 0, and is not observable if CVi = 0, where (A^, 14, Ui) 
is the eigenvalue-eigenvector pair. 

We partition the eigenvalues and eigenvectors (A, V, U) into 
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four parts; 


A = 



^ Aco 

V 



(4) 


where {Aco,Vco,Uco) are the controllable and observable 
modes, (Acs, 14g, Ucb) are the controllable but unobservable 
modes, (Ago, V"go, C/go) are the uncontrollable but observable 
modes, and (Agg, Vso, Ubb) are the uncontrollable and unob¬ 
servable modes. 

In this paper, we consider the model reduction problem for 
large-scale systems with a large number of inputs/outputs. The 
goal is to construct an ROM such that the outputs of the ROM 
Ur are close to the outputs of the full order system y, i.e., 
\y — yr \ is small. 

Denote the number of the controllable and observable modes 
is exactly I throughout the paper. First, we summarize all the 
assumptions made in this paper. 

• Al. A is stable and diagonalizable. 

• A2. Z < N. 

• A3. U',,B = 0, U',sB = 0, CV,o = 0, CVbb = 0. 

• A4. U^^B = eCi,U'^sB = eC 2 ,CVco = eC^.CVBo = 
eC 4 , where e is a small number, Ci, 02 , Ca, C 4 are con¬ 
stant matrices. And if ||t/co^ll = ^(IICsH), HCVgoU = 
OdlCsIl), then 11^111,110211,110311,110411 = 0 (|| 05 ||), 
where O 5 is a constant matrix. 


The structure of the paper is summarized in Figure |2] 
Discussion on Assumptions. For most of the practical 
applications we consider, A1 is satisfied. A2 needs to be 
satisfied for the system to have an ROM, this assumption 
is typically satisfied for a large-scale system. It should be 
noticed that from Definition 2, assumption A3 is the state¬ 
ment of controllability/observability of the different modes 
of the system. However, in practice, U^^B, CVcb are never 
exactly zero, and hence, in assumption A4, we assume that 
||17goi?|| oc e, IlCygoll oc e, where e is small. 


Algorithm 1 BPOD Algorithm 

1) Collect impulse response Xb of the primal system 
O: Use bi,i = 1, • • ■ ,p as initial conditions for O 
with Uk = 0. Take m snapshots at discrete times 
f 15 f 25 * *' 1 fm. and 

[^1 (^ 1 ); ‘ (Z'l); ‘ * * 7^1 (^m); ‘ ' i ^p (^m)] 5 (5) 

where Xi{tk) is the state snapshot at time tk with bi 
as the initial condition, k = 1,2, ■ ■ ■ ,m and i = 

1,2,- ,p. 

2) Collect impulse response Zb of the adjoint system (|2]): 

Use c', j = 1, • ■ • ,q as initial conditions for (O with 
Vk = 0- Take n snapshots at time step ,tn, 

and 

Zb = [zi{ti),--- ,zi({„),••• ,Zg{in)], (6) 

where Zj(tk) is the state snapshot of the adjoint sys¬ 
tem at time with c' as the initial condition, k = 

- ,n and j = 1 , 2 ,-- - ,g. 

3) Construct Hankel matrix Hb'. 

Hb = Z'bXb, (7) 

4) Solve the SVD problem of Hb'. 

where Eg contains the first I non-zero singular values, 
and {Lb, Rb) are the corresponding left and right singu¬ 
lar vectors. Ei contains the rest of the singular values. 

5) Construct the BPOD bases: 

Tb = Sb = (9) 

6 ) The ROM is: 

Ab = SbATb, Bb = SbB, Cb = CTg. (10) 


III. Simplified Analysis 

In this section, first, we briefly review the BPOD algorithm, 
and then illustrate in a simplified fashion how the transforma¬ 
tion bases and the Markov parameters of the ROM constructed 



Fig. 2. Structure of the Paper 


using the BPOD algorithm can be related to the controllable 
and observable modes (Ago, Vco, Uco) of the system matrix A. 
The simplified analysis is critical to understand the intuition 
behind the proposed RPOD* algorithm in Section HVl 

A. BPOD Algorithm 

Consider the stable linear system ([I]l-(|2]l, let B = 
[5i, • • ■ , bp] where bi,i = 1, • • • ,p is the column of B, 
and C = [c'l,--- ,Cg], where = 1, • • • ,qis the row 
of C. The BPOD Algorithm 0 is summarized in Algorithm 

m 

B. Simplified Analysis 

First, we construct a modal BPOD ROM by projecting 
the BPOD bases {Tb, Sb) onto the ROM eigenspace as in 
Algorithm |2] 

Under assumptions Al, A2, and A3, we have the following 
result. 
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Algorithm 2 BPOD modal ROM Algorithm 


1) 

Construct BPOD ROM (Ab, Bb,Cb) and BPOD bases 
(Tb,Sb) using BPOD Algorithm [T] 

2) 

Solve the eigenvalue problem of Ab'. 



Ab = PbAbPg^, 

(11) 


where Ab are the eigenvalues, and Pb are 
sponding eigenvectors. 

the corre- 

3) 

Construct BPOD modal bases: 



^b = Pb~^Sb,'i>b = TbPb, 

(12) 

4) 

The modal ROM is: 



Ab = ^bA'i/b, Bb = 4>bB, Cb = C-^b- 

(13) 


Proposition 1: Denote (Ab, Bb,Cb) as the modal ROM con¬ 
structed using Algorithm |2] ^'b) are BPOD modal bases. 

Then Ab = ^bA'^b = A^o, Bb = Cb = C-^b = 

CVco, where (Aco, Bco, Vco) are the controllable and observ¬ 
able modes of the system, and CbA\Bb = CA^B, i = 1, 2, • • •. 

Proof: Consider the snapshots in the primal snapshot 
ensemble ©, 

xAtk) = A^’^h, (14) 

where i = 1, • • • ,p and fc = 1, • • • ,m. Hence, 

{xi(tu), •••, Xp{tk)) = A*’^B, (15) 

and the snapshot ensemble Xb can be written as: 

Xb={A*^B A*^B ••• (16) 

Under assumptions A1 and A3, and from (IHi, we have: 

= (Ko Kb) (“I) , (17) 

where are the coefficient matrices. Substitute ([TtT i into 

([I6]l, and Xb can be written as: 

Xb = (Ko Ko) (""r) , (18) 

V^co / 

where the coefficient matrices. Similarly, 

^b = (Uco Ubo) , (19) 

where P^g, figg are some coefficient matrices. 

Hence, 

Z'bXb = {{PD'U'gg + {Plg)'U'gg){Vggalg + Vg-galg), 

= (PLyal, (20) 

where t/'„Ko = 0, = 0, C/l„Ks = 0. 

Under assumption A3, there are exactly I non-zero singular 


values in the resulting SVD problem, i.e., 

Z'bXb = Lb'ZbR'b, (21) 

where Sf, C are the I non-zero singular values and 

{Lb,Rb) are the corresponding left and right singular vectors. 
Moreover, 

Z'bAXb = {Plg)'Kgga\g. (22) 

Now, consider the BPOD ROM: 

Ab = SbATb = 

= Aco . (23) 

^ ^ ^ ^ 

Pb Pb 

We show that A^o are the eigenvalues of Ab, and Pb are the 
eigenvectors as follows: 

PbPb = 

= T.y^'^L'bLb'ZbR'bRbZ^y^'^ = I. (24) 

Also, 

PbPb = 

= o:lgAPl)'a\g)+{Pl)' = I, (25) 

where (.)’*' denotes the pseudoinverse of (.). Hence, Pb = Pg^ 
and from (|2^ . 

Acg = {Pp^Sb)A{TbPb). (26) 

$1, 

Also, we prove that ${, are biorthogonal as follows. 

= TbPbPg-^Sb = TbSb, 

K^b = Pb~^SbTbPb = Pp^Pb = /, (27) 

where {Tb,Sb) are BPOD bases, and are biorthogonal. 

Hence, 

^b = TbPb = XbRb^y^'^^y^'^L'b{Plg)' 

= (Vggal + Vg-galg)RbY.y^L'APlg)' 

= VcgPbPb + VggCe = Ko + VcoCe , (28) 

where Cq = aggRb'Zp^L'APcg)'■ Similarly, 

^b = Pb-"Sb = U^g + CrU(g, (29) 

where = a^gRb^p^L'APgg)'■ Under assumption A3, the 
modal BPOD ROM constructed using (thf,, $b) is: 

Ab = 4>bA4'b = Agg, 

Bb = $' i? = U'ggB + C^U'ggB = 

Cb = C'fb = CVgg + CVggCe = CVgg . (30) 

And the Markov parameters of the ROM are: 

CbAlBb = CVggAlgU'ggB = C A^ B . (31) 

■ 

Discussion on Proposition [T] Recall the impulse response 
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snapshot ensembles collected in the BPOD are: 



Nxl Ixpm Nxpm N xl Ixqn Nxqn 


(32) 


and the Hankel matrix is; 


Hb = Z'Xfc = e ( 33 ) 

qnxl lx pm 

There are two main parts to the computation: 

1) The primal and adjoint snapshot ensembles Xb G 

G and hence, the construction of 

Hb takes time 0{pqmnN). 

2) The computational cost to solve the SVD of Hb is 
0{imii{p^rn'^qri,prnq^n^}). 

However, Hb is only rank “Z”, where I N,pm,qn, and 
from the development of Proposition [T] we see that the modal 
BPOD ROM given in ( l30b is completely determined by the I 
controllable and observable modes and is invariant to the data 
Xb and Zb, i.e., as long as the snapshot ensembles can be 
written as: 



Nxm Nxl Ixrn Nxm 


+ (35) 

Nxn Nxl Ixn N xn 

where a *, /3* are rank I constant matrices, and a *, are 
some constant matrices of suitable dimensions, then under 
assumptions Al, A2, and A3, the following corollary holds. 

Corollary 1: Denote (A*, B*, C*) as the modal ROM con¬ 
structed using Algorithm 12] with snapshot ensembles X*,Z* 
as in (O. Then A* = Aco,B* = U'^oB,C* = CVao, where 
{Aco,Uco,Vco) are the controllable and observable modes of 
the system, and C*{A*yB* = CA^B, i = 1, 2, • • •. 

Corollary [T] can be proved by replacing Xb, Zb in Proposi¬ 
tion [T] with X*,Z*, and the proof is omitted here. We make 
the following observations. 

Observation: 1) The snapshot ensembles do not have to be 
collections of the impulse responses as in BPOD. 2) Only I 
snapshots may be enough to extract all the controllable and 
observable modes of the system. 

Bearing this observation in mind, in the next section, we 
introduce the RPOD* algorithm which generates the snapshot 
ensembles using exactly “Z” snapshots, such that the Corollary 
[D holds. 


IV. Computationally Optimal Randomized Proper 
Orthogonal Decomposition (RPOD*) 

In this section, first, we define the computationally optimal 
snapshot ensemble, and prove that the snapshot ensembles 
generated by perturbing the primal/adjoint system with white 
noise are computationally optimal snapshot ensembles. Then 
we propose the RPOD* algorithm in Section lTV-BI We discuss 
implementation issues of the algorithm in Section IIV-CI 


A. Computationally Optimal Snapshot Ensemble 


Under assumptions Al, A2 and A3, we define a computa¬ 
tionally optimal snapshot ensemble of the system as follows. 

Definition 3; A computationally optimal primal snapshot 
ensemble X* is an Z-snapshot ensemble of rank Z which can 
be written as in Eq. (l34l i. i.e., m = 1. 

A similar definition suffices for the computationally optimal 
adjoint snapshot ensemble Z*. 

Consider the system ([T]|-(|2|, under assumption Al that A is 
stable, there exists a hnite number tss, such that ||A*''''|| Ri 0. 
Under assumptions Al, A2 and A3, we have the following 
result. 

Proposition 2: Perturb the primal system ([1]) with white 
noise Uk, and collect m snapshots at time - ■ ■ , tm, where 
tm > tss, and IIA*'’® II w 0. Denote the snapshot ensemble as 
Xr = {xi X 2 • • • Xm)- If 771 = Z, where Z is the number 
of the controllable and observable modes of the system, then 
Xr is a computationally optimal snapshot ensemble. 

Proof: For the snapshots taken before tss, the state 
snapshot Xk at time k is: 

k-l 

Xk = A^Bu{k — i),k < tss- (36) 

i=0 


Suppose there is a snapshot ensemble Xf which takes tss 
snapshots at time k = \ to k = tss, then from ( |36] |, the 
snapshot ensemble Xf can be written as; 


Xf 


{Xl X2 


(u{l) 
0 
0 


xt^^) = {B AB ••• A*--iS) X 

Xi 


77 ( 2 ) * ■ • uitss 

u(l) • ■ • u{tss 

0 • ■ • u(tss 


1 ) 

2 ) 

3) 


ll(t'ss) ) 
u{tss - 1) 
'aitss 2) 


, (37) 


Vo 0 


0 77 ( 1 ) ; 

•v 

n 


where Xb is the BPOD snapshot ensemble from time k = 0 
to k = tss — 1- Denote U = (wr uj 2 ■ ■ ■ where tOi is 

the column of U. The U matrix above has columns that 
are linearly independent since it is upper triangular. 

Since Xr consists of m columns of X/, Xr can be written 
as: 


Xr = {xi ••• Xm)=Xb(uJi ••• IVm), (38) 

where (wi • ■ • utm) are the corresponding columns of U, 
and hence, Ui has full column rank. 

For the snapshot Xk taken after tss, Xk could also be written 
as: Xk = XbUJk, where ujk is a column vector whose entries 
are white noises. Therefore, tok is independent of wi, • • ■ , uim 
in (l38l l, and hence, for all the snapshots collected in Xr, Xr = 
XbCli, where Ui has full column rank. 
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Recall that from ( fTST i, Xf, can be written as: 

Xb = Vco Oi^co + ^coOicdi 
Nxllxptss Nxptss 

and has full row rank. Thus, when m = I, 


Xr = Xbni=V,o +Vcoals^u 


Nxl Ixptiiii ptsaXl 

— Ctno TV^oCtcoi 


(39) 


(40) 


where rank Uco = I- Hence, Xr is a computationally optimal 
snapshot ensemble. ■ 

Similarly, we take n = I snapshots by perturbing the adjoint 
system (|2]i with white noise Vk, and the adjoint snapshot 
ensemble can be written as: 


Zr — U^^P^+UcoPco, (41) 

Nxl Ixl 

where rank Pco = I, and Pco is a constant matrix. Thus, Zr is a 
computationally optimal snapshot ensemble. In Section IIV-CI 
we discuss about how to choose m, n and the snapshots in 
practice. 


B. RPOD* Algorithm 

The RPOD* algorithm is summarized in Algorithm |3] 

Under assumptions Al, A2 and A4, the following result 
holds. 

Proposition 3: Denote {Ar,Br,Cr) as the ROM con¬ 
structed using RPOD* following Algorithm [3 If we keep 
the first I non-zero singular values in (l44l i. then \\CrAl.Br — 
CA'‘B\\ (X 0(e), i = 1, • • •, where e is a small number defined 
in assumption A4. 

The proof of Proposition |3] uses perturbation theory ifTSl . 
lfT9l to extend the proof of the idealized Proposition [T] such 
that A4 holds instead of A3. Under assumption A4, the actual 
snapshot ensembles can be written as: 


Xr 


(Ko Vbo) 



(49) 


where 5aco = ectgo, 5aco = ectgs and e is small. Therefore, 
IjUgo^ago + VcoSasoW = 0{e) are small perturbations of the 
ideal snapshot ensemble, and we can expect the ideal result to 
be perturbed by a small amount as well. 

The formal proof is shown in Appendex 
Corollary 2: e is assumed to be a small number in assump¬ 
tion A4, and e can also be related to ct;+i as follows. 


\\CrAlBr - CA^B\\ oc 0{e) oc 0{ai+i). (50) 

The proof is shown in Appendex |B] 


C. Implementation Issues 

Here we discuss some implementation problems in the 
RPOD* algorithm. We give the insight into how to collect 


Algorithm 3 RPOD* Algorithm 

1) Perturb the primal system O with white noise Uk, 

collect m snapshots at time step - ■ ■ Am, where 

tm > tes, ||^*'’'’|| ~ 0, m > 1. Denote the snapshot 
ensemble Xr as: 

Xr = (Xi X2 ■■■ Xm) ■ (42) 

2) Perturb the adjoint system (|2]l with white noise Vk, 

collect n snapshots at time step An, where 

tn > tss, n > 1. Denote the adjoint snapshot ensemble 
Zr as: 

Zr = {Zi Z2 ■■■ Zn) . (43) 

3) Solve the SVD problem: 

z;Av = (t, "„)(§). w 

and truncate at ai, where I is the number of control¬ 
lable and observable modes present in the snapshot 
ensembles. Sr contains the first I non-zero singular 
values (Ji > a 2 >■■■,> cri > 0, {Rr,Lr) are the 
corresponding right and left singular vectors. 

4) Construct the POD bases: 

Tr = XrRrK^^^, S'r = V^Z[ . (45) 

5) Construct the ROM A, find the eigenvalues Ar and 
eigenvectors Pr of A. 

A = SrATr = PrArPr~\ (46) 

6) Construct new POD bases: 

$r = Br~^Sr, 'i’r = TrPr- (47) 

7) The ROM is: 

Ar = ^-rAT-r, Br = $rB, Cr = CT-r (48) 


the snapshot ensembles, and how to select the ROM size. 

Snapshot selection From the analysis in Section IIV-AI 
we only need to collect m = I snapshots from the primal 
simulations. However, I is not known a priori, thus, in practice, 
we start with a random guess m « N, where N is the 
dimension of the system, or we can choose m from experience. 
For instance, in a fluid system with 10® degrees of freedom, 
m is 0(10) ~ O(IO^). Similarly, we guess n, and then we 
check the rank of Z'^Xr- If Z'^Xr has full rank, i.e., rank 
{Z'rXr) = min (m, n), then it is possible that we did not take 
enough snapshots, and hence, we increase m, n, until rank 
(Z'Xr) < mm{m,n). 

For the primal simulation, we take m snapshots from 
one simulation with zero initial condition, and white noise 
perturbation u{k). We assume that the snapshots are taken at 
AT, 2AT • ■ • , mAT, WLOG. Here, AT is a small constant, 
and we require that mAT > tss, where || A*'*'’ || r:; 0. In Fig. [3 
we show one simulation result comparing the accuracy of the 
ROMs using AT = 3, 5,10, 20, 50 for the atmospheric disper- 
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Fig. 3. Snapshot Selection: Effect of AT 

sion problem introduced in Section IVI-BI Here tss = 900, and 
m = 300. The output relative error is defined in (l59l l. and we 
take the average of the output relative error for each ROM. It 
can be seen that as AT increases, each column in is well 
separated, and hence, the ROM is more accurate, while it takes 
longer time to generate the snapshots. Thus, this is a trade-off 
between the accuracy and the computational efficiency. 

ROM size selection In Proposition |3 we assume that there 
are I controllable and observable modes, and we keep exact 
I non-zero singular values. However, I is not known as a 
priori. We prove in ll20l that if we keep k non-zero singular 
values, and k > I, then undesired noise will be introduced. 
If k < I, not all the controllable and observable modes can 
be recovered. Therefore, in practice, we decide I by trial and 
error. Since rank (Z'^Xr) > I is always true, we start with 
k = rank {Z'^Xr), and check the eigenvalues of A = SAT. 
From the development in ll20l . we can see that when k > I, 
k — l eigenvalues of A are small, which are the perturbations of 
the zero eigenvalues. If k >> I, then the perturbations in the 
ROM is too large to be neglected, which results in unstable 
eigenvalues of A. Thus, we keep decreasing the value of k 
until A is stable. If there are some small eigenvalues which 
are approximately zero, we know k > I, and we decrease the 
value of k until we reach a region [1,1 +a], where a is a small 
number, such that most of the eigenvalues of A remain the 
same for different value of k (I controllable and observable 
modes with k — I perturbations of the zero eigenvalues), then 
we stop and pick the number I as the number of non-zero 
eigenvalues of A. 

V. Comparison with related Algorithms 

In this section, we compare the RPOD* algorithm with 
BPOD, random projection and BPOD output projection algo¬ 
rithm. 

A. Comparison with BPOD 

First, we summarize the differences between the BPOD and 
RPOD* algorithm. 

As we mentioned in Section HID p + q simulations are 
needed for BPOD algorithm, where p is the number of 
inputs and q is the number of outputs. Let Xf,, Zb denote 
the impulse response snapshot ensembles that need to be 
collected in BPOD algorithm from time step (I,-- - ,fss)- 
Thus, Xb e ifi^^P*^\Zb e It is expensive to store 

(p -I- q)tss snapshots, and it is expensive to solve the resulting 


SVD problem due to the large size of the problem. For RPOD* 
algorithm, only (1 primal + 1 adjoint) simulations are needed, 
and only m+n, where m, n <C tss snapshots need to be stored. 
Also, it is easy to solve the resulting SVD problem. 

Another practical problem with impulse responses snapshots 
is that the snapshots after some time are dominated by very few 
slow modes, and including these snapshots does not give much 
new information. On the other hand, the RPOD* trajectories 
are white noise forced, and all the modes are always present in 
all the snapshots due to the persistent excitation of the white 
noise terms. Hence, the RPOD* snapshots can be taken till tss, 
and be assured that all of the relevant modes will be captured. 


B. Comparison with Random Projection 

From the analysis in Section IIV-AI we see that the snapshot 
ensembles collected in RPOD* can be written as; 

Xr= Xbfll,Zr = Zbfl2, (51) 

where Xb € , Zb G j-jjg impulse response 

snapshot ensembles that need to be collected in the BPOD 
algorithm from time step (0, • • • ,tss — 1), and Ri 0. 

Hi e ^ptsaxm g ^qtssxn ^j.g matrices. We 

have: 

^ = Zl.Xr = n'^Z'bXbVti = (52) 

rank I rank I 

There is a signihcant difference between the proposed 
algorithm and a direct application of the random projection 
algorithm on BPOD. A direct application of the random 
projection would require to generate the Hankel matrix Hb 
(and Xb, Zb) first. However, in practice, the construction and 
the storage of the Hankel matrix is computationally prohibitive 
when N is large and the number of inputs/outputs is large. 
Also, the bottleneck of the random projection algorithm is the 
projection of Xb, Zb onto the Gaussian test matrices. In the 
proposed algorithm, the snapshot ensembles are constructed 
directly from the primal and adjoint simulations, and hence, 
the computational cost to generate the Hankel matrix and to 
project it onto the Gaussian test matrices is saved. 


C. Comparison with BPOD output projection 

As we mentioned in Section HID when the number of the 
outputs q is large, the computation of the BPOD adjoint 
simulations may not be tractable. To reduce the number of the 
outputs, the output projection method in ||3], ET\ is proposed. 
In this section, we compare the RPOD* algorithm with BPOD 
output projection algorithm. First, we briefly review the BPOD 
output projection method in Algorithm |4] 

In the following, we compare the BPOD output projection 
method with RPOD* algorithm. From (|5^ . the adjoint snap- 








Algorithm 4 BPOD output projection algorithm 

1) Collect the primal snapshot ensemble Xi, in (5). 

2) Compute the POD modes of the dataset Yf, = CX),. 

Y^Yb(l)k = Afc(/)fc, Ai > • • • > A„ > 0, (53) 

where Afc are the eigenvalues, and (pk are the corre¬ 
sponding eigenvectors. Thus, the POD modes 0^ = 
[^ 1 , • • • , 0s], where s is the rank of the output projec¬ 
tion, and s q. 

3) Collect the impulse responses of the adjoint system: 

Zk+i = A'Zk + C'<dsV,w = B'zk- (54) 

The adjoint snapshot ensemble is denoted as Zg. 

4) Construct the Hankel matrix: 

Hs = Z'gXb. (55) 

5) Solve the SVD problem of Hg, and construct the BPOD 
output projection bases as (Tg, Sg) using equations (8) 
and (9). The ROM is: 

As = SgATg, Bg = SgB, Cg = CTg. (56) 


shote ensemble Zg can be written as: 

Zs = (C'0s A'C"0s ■ 


{A'f“C'Qg) 


Nxsts 


= (C' A'C ••• {A'f^^C) 


Zh 


/0« 


V 


0s 


0 s 


= Zfc ^ . (57) 

N xqtss Xsiss 

And hence, the projected Hankel matrix Hs is: 

Hg = Q'Z^Xb = e'Hb. (58) 


Recall that the adjoint snapshot ensembles collected in RPOD* 
can be written as Zr = Zb^ 2 , and the projected Hankel matrix 
in RPOD* is Hr = Thus, we make the following 

remark. 

Remark 1: Both the BPOD output projection and RPOD* 
algorithms can be viewed as projecting the full order Hankel 
matrix onto a reduced order Hankel matrix with projection 
matrices 0 and ( 01 , 02 ). 

1) Differences between two algorithms: First, the infor¬ 
mation preserved in Hg, Hr are not the same. The output 
projection Pg = 0s 0s minimizes the 2-norm of the difference 
between the original transfer function and the output-projected 
transfer function, i.e. \\CA^B — QgQ'gCA'^B\\,i = I,-- - , 
is minimized. Thus, the controllable and observable modes 
preserved in Hg are an approximation of those in H},, while 
in RPOD* algorithm, the exact controllable and observable 
modes are preserved using the Gaussian random projection 
matrix. As mentioned in ED, when s < I, where I is the 


TABLE I. Computational Complexity Analysis for RPOD* 
AND BPOD Output Projection 


Construction of H 

RPOD* 

0{mnN) 

BPOD output projection 

0{pstlgN) 

Solve SVD 

0{m'^n) 

Oip^stig) 


number of non-zero Hankel singular values (number of con¬ 
trollable and observable modes), then only the first s Hankel 
singular values are the same as the full balanced truncations 
Hankel singular values. 

Another difference between output projection algorithm and 
RPOD* algorithm is that RPOD* algorithm can be used when 
both the number of the inputs and outputs are large, while 
output projection can be used when the number of the inputs 
or the outputs is large. When the number of inputs is large, 
we can construct an input projection using the adjoint snapshot 
ensemble, but when both the number of inputs and outputs are 
large, the construction of the projection matrix is not helpful. 

2 ) Comparison of computational cost: The comparison of 
the computational cost of the BPOD output projection algo¬ 
rithm with RPOD* algorithm is shown in Table U We collect 
m, n snapshots for RPOD* algorithm respectively, and tgg 
snapshots for BPOD output projection algorithm. N is the 
dimension of the system, p, q are the number of inputs and 
outputs respectively, the rank of the output projection is s, 
and without loss of generality, we assume p < q, m < n, and 
p < s. 

Now we compare the computation time to generate the 
snapshot ensembles. If we denote T as the time to propagate 
the primal/adjoint system once, then for BPOD output projec¬ 
tion algorithm, {p + s) simulations are needed, and for each 
simulation, we need to collect the snapshots up to f = tgg. 
Thus, the total computation time to generate the snapshot 
ensembles is {p + s)tggT. 

The RPOD* algorithm needs 1 primal and 1 adjoint sim¬ 
ulation till time mAT, nAT respectively, where mAT > 
tgg, nAT > tgg. When we choose AT such that mAT = 
tgg, nAT = tgg, the performance of the ROM con¬ 
structed using RPOD* is better than those constructed using 
BPOD/BPOD output projection algorithm. While, the total 
computation time to generate the snapshot ensembles is 2tggT, 
which means that RPOD* computational cost is the same as 
BPOD in a single input single output (SISO) system. The 
comparison of the computation time to generate the snapshot 
ensembles is shown in Table HIH for two examples. 


VI. Computational Results 
In this section, we show the comparison of RPOD* with 
BPOD for two examples: a one-dimensional heat transfer prob¬ 
lem, and a three dimensional atmospheric dispersion problem 
to illustrate the proposed algorithm. 

First, we dehne the output relative error: 


E, 


output 


ll^rue 

WYtrueW 


(59) 


where Ytme the outputs of the full order system, and Yrom 
are the outputs of the reduced order system. 
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The frequency response for multiple input multiple output 
systems can be represented by plotting the maximum singular 
value of the transfer function matrix max{a{H{juj))) as a 
function of frequency w. We define the frequency responses 
error as: 

Efreijuj) = I max(cr(iJt™e(ja;))) - max(cr(Tfrom(jw)))|, ( 60 ) 

where HtmeU^) is the transfer function of the full order 
system, and Hrom{j^) is the transfer function of the ROM. 

For the two experiments below, there are design parameters 
that need to be chosen manually, such as the rank of the output 
projection, when to take the snapshots, and the size of the 
ROM. We have shown the results of the best selections of these 
parameters for the BPOD/BPOD output projection algorithms 
and the RPOD* algorithm. 


A. Heat Transfer 

The equation for heat transfer by conduction along a slab is 
given by the partial differential equation: 


dT d^T ^ 

'm ~ ^ 


( 61 ) 


with boundary conditions 


m 

■< ■ 

o 



True value 
□ BPOD 

RPOD 

BPOD projection 






-RPOD 

-BPOD 

-BPOD projection 


(a) Comparison of Markov Parameters (b) Comparison of Output Relative 

Errors 


Fig. 4. Heat Transfer Problem 




(a) Compaiison of Frequency Re- (b) Comparison of Frequency Re¬ 
sponses sponses Errors 


T|a;=0 = 0, ■^|a;=L — 0, (62) 

where a is the thermal diffusivity, and / is the forcing. 

Two point sources are located at a: = 0.15m and x = 0.45m. 
The system is discretized using finite difference method, and 
there are 100 grids which are equally spaced. We take full 
field measurements, i.e., measurements at every node. The 
parameters of the system are summarized in Table |II] In the 
following, we compare the ROM constructed using RPOD*, 
BPOD and BPOD output projection. 

For RPOD*, the system is perturbed by the white noise with 
distribution N(0, 12 x 2 )- At time tss = 3000s, 0, 

thus, for the primal/adjoint simulation, one realization is 
needed, and we collect 80 equispaced snapshots during time 
t G [0, 3200]s. For BPOD, we need p = 2 realizations for 
the primal simulation, and q = 100 realizations for the adjoint 
simulations. In general, 3000 equispaced snapshots between 
[0, 3000] s should be taken in each realization for BPOD/output 
projection. However, due to the memory limits on the platform, 
only 400 equispaced snapshots can be taken and for the optimal 
performance of BPOD/output projection, the 400 equispaced 
snapshots are taken from time / G [0,400] s (first 400 dominant 
impulse responses). The rank of the output projection is 40, 
and hence, for the BPOD output projection algorithm, 40 
realizations are needed for the adjoint simulations. RPOD* 
algorithm and BPOD algorithm extract 70 modes, and BPOD 
output projection algorithm extract 50 modes. Here, we take 80 
snapshots by trial and error following the procedure in Section 

ttsa 

In Fig. St a), we compare the norm of the Markov parameters 
of the ROM constructed using three algorithms with the 
true Markov parameters of the full order system. We perturb 


Fig. 5. Heat Transfer Problem 


the system with random Gaussian noise, and compare the 
output relative errors in Fig. Hfb). In Fig. |3a), we compare 
the frequency responses of the ROM constructed using three 
algorithms with the true frequency responses. In Fig. HJb), we 
plot the error of the maximum singular value of the input- 
output model as a function of frequency. 

For all three methods, we can see that the Markov param¬ 
eters of the ROM are close to the true Markov parameters 
of the full order system. From Fig. Hfb), it can be seen that 
all three methods are accurate enough. The output relative 
error of BPOD output projection is less than 0.01%, and the 
performance of the RPOD* algorithm is much better than the 
BPOD output projection algorithm, the performance of RPOD* 
algorithm is better than BPOD algorithm because we do not 
take all the snapshots up to tss due to the memory limits. From 
Fig . H] we can see that the frequency responses error of RPOD* 
algorithm is smaller than BPOD and BPOD output projec¬ 
tion algorithm in low frequencies. With the increase of the 
frequency, the BPOD algorithm performs better than RPOD* 
algorithm, however, the errors are below 10“^°, and hence, the 
difference is negligible. The comparison of computational time 
of RPOD* and BPOD output projection algorithm is shown in 
Section rVl-CI We can see that the construction of the snapshot 
ensembles using RPOD* takes almost the same time as BPOD 
output projection, while the dominant computational cost is 
solving the SVD problem, and it can be seen that RPOD* is 
about 24 times faster than BPOD output projection. 
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TABLE II. Parameters of system 




Heat 

Atmospheric Dispersion 

System Parameters 

Domain 

X G [0, 1J(to) 

X G [0,2000J(m), y G [—100,400J(m), 
z G [0, 50] (m) 

Dimension of system 

X = 100 

O 

O 

t-H 

II 

Parameters 

a = 4.2 X 10“°(m^/s) 

Wind velocity u = (4m/s, 0,0) 

Number of Inputs 

2 

10 

Number of Outputs 

100 

810 

output projection V.S. RPOD* 

Primal Snapshots 

400 V.S. 80 

200 V.S. 400 

Adjoint Snapshots 

400 V.S. 80 

200 V.S. 400 

Snapshots taken during 

t G [0,400sJ V.S. t G 
[Os, 3200s] 

t G [0,200]s V.S. t G [0,4000]s 

ROM modes 

50 V.S. 70 

200 V.S. 380 

Hankel matrix 

16000 X 800 V.S. 80 x 80 

4000 X 2000 V.S. 400 x 400 


B. Atmospheric Dispersion Problem 

The three-dimensional advection-diffusion equation describ¬ 
ing the contaminant transport in the atmosphere is: 

dc -* 

— + V • (cu) = V ■ {K(X)Vc) + Q6{X - X^), (63) 

at 

where 

c{X,t) : mass concentration at location X = {x,y,z). 

Xs = {xstUs, Zs)' location of the point source. 

u = {ucos{a), usin{a),Q): wind velocity, a is the direction 
of the wind in the horizontal plane and the wind velocity is 
aligned with the positive T-axis when a = 0, u > 0 is constant. 

Q: contaminant emitted rate. 

V: gradient operator. 

K{X) = diag{Kx{x),Ky{x),Kz{x)) : diagonal matrix 
whose entries are the turbulent eddy diffusivities. In general 
K{X) is a function of the downwind distance x only. 

In practice, the wind velocity is sufficiently large that the 
diffusion in the ai-direction is much smaller than advection, 
and hence, assume that the term K^d^c can be neglected. 

Define a^ix) = Ky{r])dr], a^x) = K^{r])dr], 

where ay{x) = ayx{l + byx)^'^, cjz{x) = azx{l -\- bzx)^'^, 
and ay = 0.008, by = 0.00001, = 0.006, bz = 0.00015. 

The boundary conditions are: 

c(0, y, z) = 0, c(cx), y, z) = 0, c(ai, ±cx), z) = 0, 

3c 

c(a;,y,oo) = 0, —(x, y,0) = 0. (64) 

The system is discretized using finite difference method, and 
there are 100 x 100 x 10 grids which are equally spaced. The 
parameters are summarized in Table [H] There are 10 point 
sources which are shown in Fig. |6] We take the full field 
measurements (except the nodes on x = 0, oo and y = ±oo). 
In Fig. |6] we show the actual concentration of the full field at 
time t = 200s with Q as Gaussian white noise where sources 
are the dotted points in the figure. 

In this example, since the system dimension is N = 10®, 
constructing the ROM with the full field measurements us¬ 
ing BPOD is computationally impossible, and thus, we only 
compare the RPOD* algorithm with BPOD output projection 
algorithm. 


Full Order Model 



Fig. 6. Concentration Field at t = 200s 


For RPOD* algorithm, we collect the snapshots sequentially. 
The system is perturbed by white noise with distribution 
A^(0,/ioxio)- One primal simulation and one adjoint simu¬ 
lation are needed. We collect 400 equispaced snapshots from 
time t € [0,4000]s, where at time = 4000, ||^*''“|| ~ 0, 
and extract 380 modes. For BPOD output projection algorithm, 
we collect the impulse responses from the primal simulations, 
and p = 10 realizations are needed. Same as the heat 
example, we could not collect all the impulse responses from 
t S [0,4000]s due to the memory limits on the platform. 
Hence, for the best performance available in this example, 
we collect 200 equispaced snapshots from t G [0,200]s for 
each primal simulation. The rank of the output projection is 
80, and hence, 80 adjoint simulations are needed, and in the 
adjoint simulations, we collect 50 equispaced snapshots from 
t G [0, 50]s. We can extract 200 modes. The parameters are 
summarized in Table In] 

In Fig.|2ta), we compare the Markov parameters of the ROM 
constructed using RPOD* and BPOD output projection with 
the full order system. Also, we perturb the system with random 
Gaussian noise, and compare the output relative errors in Fig. 
|7Jb). The comparison of the computational time is shown in 
Section IVI-CI 

In Fig. da), we compare the frequency responses of the 
ROM constructed using RPOD* and BPOD output projection 
with the full order system. We can see that the frequency 
responses of the ROMs are almost the same as the frequency 
responses of the full order system. In Fig. db), we show the 
errors between the frequency responses. 

The comparison of the computational time is shown in Sec- 
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0 

D True value 

BPOD projection 

0 RPOD 

a 




Q 




(a) Comparison of Markov Parameters (b) Comparison of Output Relative 

Errors 


Fig. 7. Atmospheric Dispersion Problem 




(a) Comparison of Frequency Re- (b) Compaiison of Frequency Re¬ 
sponses sponses Errors 

Fig. 8. Atmospheric Dispersion Problem 


tion IVI-Cl It can be seen that for this example, the construction 
of the snapshot ensembles using RPOD* is faster than the 
BPOD output projection, and the dominant computation cost 
is the construction of Z'X, where RPOD* algorithm is almost 
16500 times faster than BPOD output projection. 

From the examples above, we can see that for both examples 
showed in this paper, the RPOD* algorithm is much faster than 
BPOD/BPOD output projection algorithm, and is much more 
accurate than BPOD output projection algorithm. 


C. Comparison of Computational Time 

Comparison of computational time is shown in Table [111] 
for two examples. All of the experiments reported in this 
paper were performed using Matlab 2013b on a Dell OptiPlex 
9020, Intel(R) Core (TM) i7-4770CPU, 3.40GHz, 4GB RAM 
machine. 


VII. Conclusion 

In this paper, we have introduced a computationally optimal 
randomized POD procedure for the extraction of ROMs for 
large scale systems such as those governed by PDFs. The ROM 
is constructed by perturbing the primal and adjoint system 
with Gaussian white noise, where the computational cost to 
construct the snapshot ensembles is saved when compared to 
perturbing the primal and adjoint system with impulses in 
BPOD/BPOD output projection algorithm. Also, it leads to 
a much smaller SVD problem, and an orders of magnitude 
reduction in the computation required for constructing ROMs 
over the BPOD/ BPOD output projection procedure. The 
computational results show that the accuracy of the RPOD* 


is much more accurate than the BPOD output projection 
algorithm. 


Appendix A 

Proof of Proposition [3] 

Proof: The adjoint snapshot ensemble Zr can be written 
as: 


Zr — UcoPco + UcoSPco + UcoPco + UcoSPco, (65) 

where 6f3co = f-Pco, ^/3eo = e/3eo, e is dehned in assumption 
A4, PcojPco are coefficient matrices. From ( |49] | and (l65T l. 

Zj,Xr — Pco^co “b “b Pco^^co “b ^j^co^^co 

— P'co^co + r {P'co^co + P'co^co) +0{e'^), 

" -V-' 

El 

— P'co^co + iEi + O(e^). (66) 


And similarly, 

Z^AXr — j^roAcoGCco “b 6 {Pco-^coC^cd ~b PqqAcoO^co) “bO(e ) 

V ^ 

E 2 


= P'^^kroaroEeE2 + 0{t^). (67) 

There are controllable and observable modes, then rank 
{Z'^Xr) > I due to the small perturbations. Denote 

Hr = fl'roC^co — Lr^rR'r “b Lq^oRoj 
Hr = Z'^Xr = P'^gOtco + eA'l = LrZirR'r + LoEiqRo, (68) 

where fTr,Hr G and WLOG, n < m. 

Here, fir is the ideal Hankel matrix constructed with the 
simplifying assumption (assumption A3 is satished), and it 
can be seen that the true Hankel matrix Hr (assumption A4 
is satisfied) can be viewed as adding a small perturbation of 
Hr, i-S-, Hr = Hr + 

_Er_G contains the I non-zero singular values of fir and 
{Lr,Rr) are the corresponding left and right singular vectors. 
Sp G_ 5R(*‘-0x(n-i) = 0 are the rest singular values, and 
(Lo, Ro) are the corresponding left and right singular vectors. 
Similarly, E,. G contains the first I non-zeros singular 
values of Hr, and Eq G 3 fi(”“ 6 x(r !,-0 contains the rest singular 
values. The left and right singular vectors are partitioned in the 
same way. 

1 ). From the perturbation theory ifTSl . lfT9l . the perturbed 
singular values and singular vectiors {'Er, Lr_, RrJ are related 
to the singular values and singular vectors {Er, Lr, Rr) as: 

Er = Er + ^E^ -f 

Lr — Lr “b kLr, Rr — Rr “b kRr, (69) 

where E^,XLr,XRr are some matrices, and ||ALr||, 
II Ai?j.|| oc 0(e). The expression of E3 is given by the follows 

ESI. 

For (fi G Er, (strictly positive singular values) with multi¬ 
plicity t, Lt,Rt are the corresponding left and right singular 
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TABLE III. Comparison of Computational Time 



Heat 

Atmospheric Dispersion 


RPOD* 

output projection 

RPOD* 

output projection 

Generate X 

0.0148s 

0.0518s 

55.59s 

30.461s 

Generate Z 

0.0340s 

0.0864s + 0.1582s(projection) 

56.23s 

421.99s + 9.287s(projection) 

Construct Z'X 

0.0292s 

0.0461s 

0.321s 

5311s 

Solve SVD 

0.1052s 

2.5550s 

0.4859s 

9.118s 

Total time 

0.1832s 

2.8975s 

112.6269s 

5781.856s 


vectors, the perturbed singular values ai € and 

X,iR'tE[Lt + L',E,Rt) , 


cTi = ai + e- 


■ +0(e ),i — 1, 




where Ai(.) denotes the eigenvalue of (.). Thus, 


E 3 — 


V 



(71) 


where ei, • • • , e; are the coefficients computed from dTOl l. 

2). The ROM A = SrATr can be proved to be a perturbation 
of Ab in d^ . 

The POD bases T^, Sr are; 


be proved that Aco are the eigenvalues of A, and Pr are the 
corresponding eigenvectors. 

3). From perturbation theory ||22l, since 

A = A + A 3 + 0{e'^)=PrArP-\ (80) 

it can be proved that ||Pr —^r|| < ll^all = ||Ar —Aco|| < 
II A 3 II = fcae. Thus, 

= TrPr = XrRr^r^^^Pr + A 4 , 
q>r = p-^Sr = + As, (81) 

where A 4 ,A 5 are some matrices and ||A 4 II, HAsH oc 0 (e). 
Substitute dSTl i into the ROM Markov parameters, and collect 
the small perturbation terms, we have: 

CrAlBr = C'i>rAi.<^>rB = CV^oKoU'coB' + A, (82) 


Tr = XrRrK^^‘^,Sr = L'. 


(72) 


where A is some matrix, and ||A|| oc 0(e). 


Therefore, 

i = SrATr = S-P^L'^Z^AXrRrY-P^ 

= Y-P^L'AP'ro^coaco + eE2)RrY-P^. (73) 

From (TTOl l, it can be proved that; 

+ (74) 

where Cr is a diagonal coefficient matrix. Therefore, 

= (E-1/2 + tCr){K + at;) = + Ai, (75) 

where Ai is some matrix, and ||Ai ||2 = fcie, fci is a constant. 
Similarly, 

= RrY-^/^ + A 2 , (76) 

where A 2 is some matrix, and IIA 2 II 2 = ^ 26 , ^2 is a constant. 
Thus, 

A = Y-P^L'Xo^coacoRrYr'^^ + A 3 + O(e^), (77) 


Appendix B 

Prooe oe Corollary | 2 ] 

For ai € Eo (zero singular values) m with multiplicity 
n — I, the corresponding left and right singular vectors are 
Lo, Ro- The perturbed singular values ai € Eq and 

CTi = e^JXi{R'^E[LoL'^EiRA, i = 1, • • ■ n - T (83) 
From (iTOl i and (l83t . we can see that: 

cr; = CT; + e;e + 0(e^),cri+l = e;+ie, (84) 

Hence,we have: 

ai+i oc 0(e), (85) 

and 

WCrAlBr - CA^B\\ oc 0(e) oc 0{ai+i). (86) 


where A3 is some matrix, and IIA3II2 = k ^ e , is a constant. 

If we let Reeerences 


then 


A = Y-^/^L'Xo^coaroRrK^'^. 

Pr Pr ^ 

A = A -\- A3 + O(e^) = 


(78) 


(79) 


Following the same proof in Section IIII-BI ( (l2^ - (l25T l). it can 
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